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Abstract. - It is shown that the Shan-Chen (SC) model for non-ideal lattice fluids can be made 
compliant with a pseudo free-energy principle by simple addition of a gradient force, whose expres- 
sion is uniquely specified in terms of the fluid density. This additional term is numerically shown 
to provide fairly negligible effects on the system evolution during phase-separation. To the best of 
our knowledge, these important properties of the SC model were not noted before. The approach 
developed in the present work is based on a continuum analysis: further extensions, more in line 
with a discrete lattice theory (X. Shan, Phys Rev E, 77 066702 (2008)) can be envisaged for the 
future. 



Multiphase fluid flows arc ubiquitous in nature, sci- 
ence and technology [1,2]. Owing to the complexity of 
the underlying physics, often leading to the formation of 
Complex space-time changing interfaces between the vari- 
ous phases/components, analytical methods have a limited 
range of applicability for the quantitative study of multi- 
phase flows. Consequently, flexible and efflcient methods 
for the numerical simulation of multiphase flows are in 
constant demand. In the last decade, a new class of meth- 
ods, based on a suitably discretized version of the Boltz- 
rnann kinetic equation, i.e. the lattice Boltzmann (LB) 
models [3-5] have been developed [6-14]. 

The main appeal of the lattice kinetic approach is the 
simplicity of its underlying dynamics, which consists of 
a free-streaming of a set pseudo-particles along a set of 
prescribed directions (discrete speeds), followed by a col- 
lisional relaxation around local equilibria. Crucial to the 
extension of the LB scheme to non-ideal fluids, is the possi- 
bility of encoding potential energy interactions at the level 
of a simple forcing term, acting upon the discrete Boltz- 
mann distributions. Suitable choices of this forcing term 
have proven capable of triggering a fairly rich non-ideal 
fluid behaviour, such as phase-transitions and interface 
formation/propagation [6,15,16], with no need of track- 
ing the shape of the interface, but rather having it emerge 
from the underlying mesoscopic interactions. 



A widely used lattice Boltzmann algorithm [3] with in- 
terparticlc interactions is the pseudo-potential scheme due 
to Shan-Chen (hereafter SC) model [6,7]. The evolution 
over a unit time lapse reads as follows: 

fi{x + cut+l)- fi{x,t) - (^fi{x,t) - f^'''\p,u')) 

(1) 

where fi(x,t) is the kinetic probability density function 
associated with the particle velocity c/ , r is the mean colli- 
sion time, Z/'^^'' (p, u) the local equilibrium, corresponding 
to a Maxwellian distribution with local density p and local 
fluid speed u. Finally, pu' — pu + tF where F represents 
a general forcing term encoding inter-molecular interac- 
tions. To be noted that this force enters the scheme in the 
form of a shift, u' — u = Ft/ p, of the local flow speed. 
In the original SC model [6,7], the bulk interparticle in- 
teraction is proportional to a free parameter (the ratio of 
potential to thermal energy), Q, entering the equation for 
the momentum balance: 



F^'^ix) = g^ix) > ] wi{\ci\-')i,{x + ci)c\ 



I 



(2) 



^) are static weights normalized to unity [23], 
= ip{p{x,t)) is the (pseudo) potential de- 



where ( | c; 
and ip{x,t) 

scribing the fluid-fluid interactions triggered by inhomo- 
geneities of the density proflle. 
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This pseudo-potential may also be viewed as a gener- 
alized density, obeying the general properties of going to 
zero in the limit p ^ 0, and saturating to a constant 
value at large densities. The mapping p ^(p) may 
be heuristically interpreted as the result of a rcsumma- 
tion of many-body collisions at all orders. On practical 
grounds, its main effect is to set the intcrmolecular forces 
to zero at high density, a property which wc may pic- 
torially call 'asymptotic freedom', since it indicates that 
at high-dcnsity/short-distances, particles are set free from 
mutual interactions. This property is crucial to prevent 
density collapse due to attractive interactions, and surro- 
gates the effect of hard-core repulsion, which would other- 
wise jeopardize the numerical stability of the scheme. In 
fact, the replacement of physical density with a general- 
ized is probably to be regarded as a general constraint of 
lattice models of non- ideal hydrodynamics. 

The presence of the pseudo-potential turns the bulk 
ideal equation of state Pb{p) = pT (with T the local tem- 
perature ) into a non-ideal one, Pb{p) = pT — ^tp'^- It 
is readily seen that a sufficiently large value of the cou- 
pling parameter Q leads to a Van der Waals-like loop, 
corresponding to the formation of a stable liquid-gas in- 
terface (phase-separation) [6,7]. The above concerns the 
behaviour of the bulk flow, sufficiently away from the in- 
terface. In order to discuss interface physics, one needs 
to inspect the contribution of the interactions to the non- 
ideal pressure tensor, defined by the relation: 



sc 



d,{pT). 



(3) 



The properties of this tensor have been investigated in a 
number of recent pubUcations [17,18]. It should be noted 
that a direct transcription of the force balance on a lattice 
with a finite set of velocities leads to a discrete form of 
the pressure tensor [18] that differs slightly from the con- 
tinuum definition (3). However, by suitably enlarging the 
set of discrete velocities [17], is well approximated ^ 
by its continuum limit 



P — 



n(p)-||vv^|2-|v'A^ 



5,j + ^d,^dj^p. (4) 



that satisfies (4) upon Taylor expansion of the forcing 
(2). In the sequel, we shall refer to this expression of the 
pressure tensor, leaving the discussion on the discrete 
lattice effects to a future work. 

At variance with other LB formulations for non-ideal 
fluids, the SC was not derived from a free-energy func- 
tional, which has often been pointed out as a theoretical 
liability of the model. 

Yet, in view of its practical success, it is tempting to 
argue that there might be an 'hidden' pseudo-free energy 
(PFE) behind it. In this work, we shall show that the no- 
tion of generalized density can indeed be derived from a 



suitably generalized (pseudo) free-energy functional. This 
is a potentially important result in two respects; first, it 
puts the SC method on a solid foundational basis within 
the general framework of non-equilibrium statistical me- 
chanics. Second, it opens new perspectives for the for- 
mulation of a thcrmo-hydrodynamical lattice formulation 
of non-ideal fluids. A very convenient aspect of the SC 
model over direct free-energy formulations [8] is that local 
equilibria keep a MaxwcUian shape, only shifted by the 
amount Ft/ p. This offers better control on its numerical 
stability, and also a direct link with self-consistent thermal 
equihbria [19]. 

The former property is key for lattice implementations, 
since it permits to preserve the same "stream-collide" 
structure of LB schemes for ideal gases, the effects of 
non-ideal interactions being entirely incorporated within 
the self-consistent velocity shift. It is worth reminding 
that the stream-collide structure is essential to the effi- 
ciency/accuracy of the method, since the lattice streaming 
operator (Ihs of equation (1)) is exact, and collisions are 
conservative up to machine round-off. 

Free-Energy derivation. — We begin by consider- 
ing a free-energy functional in the standard weak-gradient 
approximation [1]: 



C{p,Vp) = f{p) + -\Vp\'' 



(5) 



where the bulk free energy f{p) is connected to the bulk 
pressure via the familiar equilibrium Legendre's trans- 
form: 

Pb{p)^pf^-f{p)- (6) 

The term f |Vpp describes the energy surplus required 
to build up and maintain the interface, and is clearly re- 
lated to surface tension [1]. Despite its manifest link to 
underlying molecular interactions [20], to the best of our 
knowledge, no first-principle derivation of such relation is 
known. As a result, this term leaves us with some degree 
of latitude to design a free-energy functional consistent 
with the equilibrium states of the SC model. 

Since the off-diagonal terms of the SC pressure tensor 
in (4) take the form « di')pdjip, it is natural to postulate 
the following hybrid (p-ijj) form of free-energy functional: 



/:(p,vv) = /(p) + -|vV'l' = /(/>) + 



""f^) \yp\\ 



2 \dp J 



^Technically speaking, this limit corresponds to the case 64 >> 1 
in equations (22a) and (22b) and e fa 1 in equation (25) of [18] . 



It is worth noting that this 'hybrid' form is equivalent 
to the standard one (5), only with a density-dependent 

surface coupling, . = nip) = . (f This property has 
been recently exploited to produce a number of non-trivial 
dynamical effects, such as structural arrest and ageing, in 
the framework of binary mixtures [21]. Having introduced 
the hybrid free-energy, we next proceed to identify the re- 
sulting conserved currents through a standard application 
of Noether's theorem. 



Free Energy formulation for a class of multiphase models 



Invariance of the free-energy with respect to spatial 
translations, leads to the following conserved tensorial cur- 
rent [22] 

(7) 



B J 



with 



J^■. 



dC 



djp. 



(8) 



The next step is to explore the properties of this conserved 
current under the constraint of mass conservation and sta- 
tionarity of the free energy. To this purpose, we introduce 
a Lagrange multiplier. A, associated with mass conserva- 
tion, and re-dcfine our PFE accordingly: 



K f dijj 
^^P^ + 2 U 



(9) 



Minimization of the constrained free-energy leads to the 
following Euler-Lagrangc equations: 



- a. 



-A = 



dp \d{dip) 

from which, the value of A is obtained: 

, df{p) dij 

A = — ; K—Ayj. 

dp dp 



(10) 



(11) 



With this value of A and the use of relation (6), the con- 
served current, hereafter to be identified with the free- 
energy consistent pressure tensor, , becomes 



rX ^ pFE ^ 



K ._ 2 dtp 



Pbip) - - ^P^AiJ 



(12) 

where we observe that, as anticipated, in order to re- 
produce the off-diagonal term nditpdjip^ a contribution 
■||VV'P in the free energy is actually required. This im- 
mediately fixes K = ^. 

Connection to Shan-Chen Model. The above 
derivation bears close links to the pressure tensor associ- 
ated with the SC forcing, in its simplest instance: 

Ff^'ix) ^ g^{x) wi{\ci\^)^{x + ci)c]. 
I 

In particular, on the assumption of gentle interfaces, hence 
neglecting higher order terms, this force (we remind that 
in our notation G > 0) reads as: 



SC 



Qd, 



1 



(13) 



where we have made use of the usual lattice tensor prop- 
erties (see [17] for a detailed description). The first term 
in (13) is related to the bulk pressure via the relation (3), 



Pb{p)^pT--gi^{p). 



The required function /(p) is readily derived from its def- 
inition as a Legendre's transform of the bulk pressure 



.f(p) 



"Ml 



di 



Pa 



where the constant of integration po affects the total free 
energy through an additive constant, with no effect on 
the time derivative, because mass is assumed to be con- 
served. The surface term ^QipdiAip has to be handled 
with care. The forcing term associated with the consis- 
tent free-energy pressure tensor in the continuum, reads 
(after setting k = Q/2 in (4)) as follows: 



P 



FE 



Pbip) 



IW-'I 



-diipdjip. 



4' 2^dp ^\ z 

(14) 

This does not match the SC formulation given in (4), 
due to the presence of the term ^p^Aip. The reason 
is clear: physical mass, associated with p, is conserved, 
while generalized mass, associated with ■0, is not. This is 
why, strictly speaking, the SC model cannot be derived 
from a free-energy functional. Nevertheless, the interest- 
ing observation is that the only difference is confined to 
the diagonal component of the pressure tensor, while the 
off-diagonal term, ^dttl^djip, is exactly reproduced. This 
is a very welcome property, since it allows to match the 
two descriptions by simply adding an ad-hoc correction 
in the form of the gradient of a -known- potential. More 
precisely, we start from the following relation: 



pFE ^ pSC , p(add) 



where we have defined: 



In the above. 



C(P) 



[C(P)A^] 4 



dp 



(15) 



(16) 



(17) 



is the Legendre's transform of the generalized versus 
physical density. This term is responsible for the non- 
conservation of the generalized mass. Upon differentiation 
of expression (15) and use of (3), we obtain 

-9,/^^^ + d,{p) = Ff^ - d.V^^{p, A^) (18) 

where we have defined 



F^^(p,A^) 



[C(P)AV'] 



(19) 



Clearly F^-^(p, A-0) is the precise free-energy potential 
needed to recover a consistent free-energy treatment. In 
other words, 

pFE ^ pSC _ Q^yFE (20) 

is the right forcing, stemming from a free-energy consistent 
equilibrium. In conclusion, the configuration minimizing 
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the free-energy functional, whose density is 



Pa 



(21) 



and the static equilibrium associated with the forcing: 

Ff^{x) = gi;ix) J2 wiiWil^Mx + ci)cl - (22) 
I 

share the same physics. 

Continuum MsLXwell Rule. Next, we proceed to 
show that the integral constraints associated with the con- 
tinuum free-energy density, lead to the Maxwell equal-area 
rule fixing the liquid and gas density (pi and pg), indepen- 
dently of the structure of the interface, i.e. independently 
on whether the physical or generalized densities are used. 
This is again crucial to the practical success of the SC 
model. To this end, we begin by minimizing C^{p,'^ip) 
in (9) along a flat and stable liquid-gas interface, where 
inhomogeneities are confined along the y coordinate. This 
sets the value of A 



A = 



dfjp) 
dp 



(23) 



with the bulk phases at the same chemical potential (be- 
dp ^yy ^ 



cause ^dyy-ijj = in the bulk) 



dfip) 



dp 



dfip) 



dp 



Upon multiplying by ^ and integrating between y 
(conventionally, the bulk gas phase) and a generic y, at 
density p{y), we obtain: 



(24) 







dl 
dy 



dif} d'^tp 
dy dy"^ 



A— 
dy' 



(25) 



Integrating once more along y, and making use of the re- 
lation (23), leads to the following identity 



Pbip)-Pb{Pg 



2^ dp 



dipy 1 
dy J p 



(26) 



By calhng Pb{pg) = Pq the constant value we wish to 
achieve, we finally obtain 



do 

iP,{p)-Po)4 



K / dip 
2 [di^ 



1 

P' 



(27) 



Since the rhs is zero in the bulk phase, we finally conclude 
that 



r\p,[p)^p,)%=Q 



(28) 



which is the desired expression of Maxwell's equal-area 
law, q.e.d. 



Numerical tests. — In this section, our theoretical 
findings are tested against numerical simulations. For the 
sake of simplicity, we shall refer to phase-separation in a 
2d system at a constant temperature (T = 1), as explained 
in [23]. 

We also have chosen^, the pseudo-potential in the form 
'\l)[p) — , which sets the free-energy density to 



fip) = pTlogp - ^pV'^ 



(29) 



where the first term is the ideal gas contribution, supple- 
mented by an interaction term. The value of Q has been 
chosen equal to Q = 4.3, so as to produce stable interfaces 
with densities pg = 0.26 (gas phase) and pi = 1.03 (hquid 
phase) approximately with an interface of the order of ten 
grid points. 

In figure 1 we show the time evolution of the total 
pseudo-free energy F{t) ~ J C{p, Vip)dxdy, under the 
SC dynamics. We have used periodic boundary condi- 
tions and an initial state with small random perturbation 
around an average density value (see plot (a) of this fig- 
ure). The PFE shows a clear monotonic decrease in time, 
indicating that the SC dynamics does indeed support some 
minimum principle. This finding is further corroborated 
by a direct comparison of the SC dynamics, as given by 
equation (2), against the PFE dynamics, as obtained from 
the corrected model (22). The numerical data clear indi- 
cate that the lack of the correction term, not only does not 
spoil the monotonic decrease of the PFE, but also provides 
a quantitative agreement with the correct PFE dynamics. 
We conclude that, notwithstanding its conceptual impor- 
tance, the term restoring compliance of the SC model with 
a continuum free-energy functional, is fairly negligible on 
practical grounds. This is the basic reason behind the suc- 
cess of the SC model, and one that clearly hinges on the 
fact that the departure between physical and generalized 
densities, as measured by the 'metric' C(p)j remains suf- 
ficiently "small" all along the evolution. The picture is 
further clarified by monitoring the detailed dynamics of 
an 'elementary' coalescence event, i.e. two droplets merg- 
ing into a single one (equilibrium state). From the results 
reported in figure 2, the trend to a minimum PFE is again 
clearly visible, if only with a much smaller change of the 
ratio F{t)/F{0), as compared to the phase-separation test, 
which was clearly involving several merging events. In par- 
ticular, the figure reveals that virtually all of the change in 
PFE is achieved during the merging event. Again, close 
agreement between the SC versus PFE dynamics is ob- 
served, indicating that the quantitative match observed 
in the phase-separation test applies to each merger event 
separately, and is not due to some form of statistical can- 
cellation between different merging events. 

Conclusions and outlook. — Summarizing, we have 
presented a number of theoretical findings on the proper- 



^with this choice, the integral for /(p) in (21) can be computed 
exactly 
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(c) t=500 (d) t=750 



Fig. 1: Numerical simulation of two-dimensional phase-separation using the SC model. Separation proceeds from 
an initial (see plot (a)) random configuration. We monitor the behaviour of the total PFE F{t) = f C{p, \/ip)dxdy 
with £ given in equations (29) and (21). The total PFE is normalized to its initial absolute value. A sequence of 
density contour plots are displayed in (a)-(d), corresponding to t = 0, 250, 500, 750 time units. In the left panel, 
we compare two characteristic trajectories of the total PFE, as obtained with a PFE consistent forcing (22) and 
with a bare SC forcing (2), with Q = 4.3 and 4^{p) — . In both cases we have used the lattice Boltzmann 
equation (1) with r = 0.7. 
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Fig. 2: Numerical simulation of two-dimensional droplet merging with the SC model. We monitor the behaviour 
of the total PFE, as computed from F{t) = / J~-{p, Vil>)dxdy with C given in equations (29) and (21). As in figure 
1, we compare two characteristic behaviours obtained with a PFE-consistent forcing (22) and with a bare SC 
forcing (2), both with Q = 4.3 and a pseudopotential given by V'(p) = ■ In both cases we have used the 
lattice Boltzmann equation (1) with r = 0.7. 



ties of the SC model for non-ideal fluids. In particular, we 
have shown that the SC model supports a genuine pseudo- 
free-energy, monotonically decreasing in time as the fluid 
phase-separates. According to the common tenet, the SC 
dynamics does not derive from a continuum free-energy 
functional; however, we have shown that it can be made 
compliant with a free-energy functional, by simply adding 
a gradient force, whose expression is uniquely specified 
in terms of the density field. Numerical simulations of a 
phase-separating fluid, provide a clear indication that the 
additional gradient force plays a fairly negligible role on 
the evolution of the system. This unveils one of the main 



reasons behind the practical success of the SC method; 
although it was not derived from a free energy minimiza- 
tion principle, the SC dynamics docs possess a pseudo-free 
energy nonetheless. To the best of our knowledge, these 
crucial properties of the SC model were not noticed before. 

The list of the good news completed, a few words of cau- 
tion are in order. First, we have related the SC dynamics 
to a continuum free energy functional. While this is in or- 
der as a zero-th approximation, the 'discrete' limit of this 
analysis needs to be inspected very carefully [18]. In par- 
ticular, one would like to relate the SC dynamics to a dis- 
crete free-energy, i.e. a free-energy functional constructed 
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directly on the lattice. The inclusion of (lattice) free en- 
ergies and, more generally, interaction energies, faces with 
the important issue of the existence of a H-theorem for lat- 
tice models with non-ideal interactions, as well as with the 
correct formulation of lattice thermo-hydrodynamic mod- 
els [19]. The investigation of lattice models with non- ideal 
interactions and thermal dynamics, represents a major is- 
sue for future research in the field, with a potential far- 
reaching impact for non-equilibrium statistical mechan- 
ics [24] (flowing systems far from equilibrium) . 

* * * 

We gratefully acknowledge discussions with R. Bcnzi, L. 
Biferale and G. De Frisco. 
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